% Show variations in cases
[~, sidx] = sort(countries);
ignore_idx = [51, 52, 53, 55, 56];
if sum(ignore_idx == cid)>0
x_dat = datetime(2020, 1, 23)+ caldays(2:size(data_4, 2));
p_b = plot(x_dat, diff(data_4(cid, :)), 'r'); hold on;
p_b.Color = 1 - (1-p_b.Color)*0.3;
plot(x_dat, diff(data_4_s(cid, :)), 'r', 'LineWidth', 2); hold on;
xx = diff(squeeze(net_infec_A(:, cid, :))')';
x_dat = datetime(2020, 1, 23)+ caldays(size(data_4, 2)+1:(size(data_4, 2)+size(xx, 2)));
X_plot = [x_dat, fliplr(x_dat)];
Y_plot = [min(xx), fliplr(max(xx))];
plot(x_dat, mean(xx, 1), 'g', 'LineWidth', 2); hold on
fill(X_plot, Y_plot , 1,....
xx = diff(squeeze(net_infec_B(:, cid, :))')';
x_dat = datetime(2020, 1, 23)+ caldays(size(data_4, 2)+1:(size(data_4, 2)+size(xx, 2)));
X_plot = [x_dat, fliplr(x_dat)];
Y_plot = [min(xx), fliplr(max(xx))];
plot(x_dat, mean(xx, 1), 'b', 'LineWidth', 2); hold on
fill(X_plot, Y_plot , 1,....
xx = diff(squeeze(net_infec_C(:, cid, :))')';
x_dat = datetime(2020, 1, 23)+ caldays(size(data_4, 2)+1:(size(data_4, 2)+size(xx, 2)));
X_plot = [x_dat, fliplr(x_dat)];
Y_plot = [min(xx), fliplr(max(xx))];
plot(x_dat, mean(xx, 1), 'c', 'LineWidth', 2, 'LineStyle', '--'); hold on
fill(X_plot, Y_plot , 1,....
% xx = diff(squeeze(net_infec_D(:, cid, :))')';
% x_dat = datetime(2020, 1, 23)+ caldays(size(data_4, 2)+1:(size(data_4, 2)+size(xx, 2)));
% X_plot = [x_dat, fliplr(x_dat)];
% Y_plot = [min(xx), fliplr(max(xx))];
% plot(x_dat, mean(xx, 1), 'm', 'LineWidth', 2, 'LineStyle', '--'); hold on
% fill(X_plot, Y_plot , 1,....
% 'edgecolor','none', ...
h(1) = plot(NaN,'r', 'LineWidth', 2);
h(2) = plot(NaN,'g', 'LineWidth', 2);
h(3) = plot(NaN,'b', 'LineWidth', 2);
h(4) = plot(NaN, 'c', 'LineWidth', 2);
h(5) = plot(NaN, 'm', 'LineWidth', 2, 'LineStyle', ':');
'Scenario 1: No strain, high vaccine efficacy', 'Scenario 2: High Vaccine Efficacy + Strain', 'Scenario 3: Fatigue/Hesitancy+No Strain', ...
'Scenario 4: Fatigue/Hesitancy+ Strain', 'location', 'northwest');
ylabel('New Reported Cases');
set(gca, 'FontSize', 18);
set(gcf,'position',[10,10,1200,450]);
x_dat = datetime(2020, 1, 23)+ caldays(2:size(data_4, 2));
p_b = plot(x_dat, diff(deaths(cid, :)), 'r'); hold on;
p_b.Color = 1 - (1-p_b.Color)*0.3;
plot(x_dat, diff(deaths_s(cid, :)), 'r', 'LineWidth', 2); hold on;
xx = diff(squeeze(net_death_A(:, cid, :))')';
x_dat = datetime(2020, 1, 23)+ caldays(size(data_4, 2)+1:(size(data_4, 2)+size(xx, 2)));
X_plot = [x_dat, fliplr(x_dat)];
Y_plot = [min(xx), fliplr(max(xx))];
plot(x_dat, mean(xx, 1), 'g', 'LineWidth', 2); hold on
fill(X_plot, Y_plot , 1,....
xx = diff(squeeze(net_death_B(:, cid, :))')';
x_dat = datetime(2020, 1, 23)+ caldays(size(data_4, 2)+1:(size(data_4, 2)+size(xx, 2)));
X_plot = [x_dat, fliplr(x_dat)];
Y_plot = [min(xx), fliplr(max(xx))];
plot(x_dat, mean(xx, 1), 'b', 'LineWidth', 2); hold on
fill(X_plot, Y_plot , 1,....
xx = diff(squeeze(net_death_C(:, cid, :))')';
x_dat = datetime(2020, 1, 23)+ caldays(size(data_4, 2)+1:(size(data_4, 2)+size(xx, 2)));
X_plot = [x_dat, fliplr(x_dat)];
Y_plot = [min(xx), fliplr(max(xx))];
plot(x_dat, mean(xx, 1), 'c', 'LineWidth', 2, 'LineStyle', '--'); hold on
fill(X_plot, Y_plot , 1,....
% xx = diff(squeeze(net_death_D(:, cid, :))')';
% x_dat = datetime(2020, 1, 23)+ caldays(size(data_4, 2)+1:(size(data_4, 2)+size(xx, 2)));
% X_plot = [x_dat, fliplr(x_dat)];
% Y_plot = [min(xx), fliplr(max(xx))];
% plot(x_dat, mean(xx, 1), 'm', 'LineWidth', 2, 'LineStyle', '--'); hold on
% fill(X_plot, Y_plot , 1,....
% 'edgecolor','none', ...
h(1) = plot(NaN,'r', 'LineWidth', 2);
h(2) = plot(NaN,'g', 'LineWidth', 2);
h(3) = plot(NaN,'b', 'LineWidth', 2);
h(4) = plot(NaN, 'c', 'LineWidth', 2);
h(5) = plot(NaN, 'm', 'LineWidth', 2, 'LineStyle', ':');
'Scenario 1: No strain, high vaccine efficacy', 'Scenario 2: High Vaccine Efficacy + Strain', 'Scenario 3: Fatigue/Hesitancy+No Strain', ...
'Scenario 4: Fatigue/Hesitancy+ Strain', 'location', 'northwest');
ylabel('New Reported Deaths');
set(gca, 'FontSize', 18);
set(gcf,'position',[10,10,1200,450]);